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FOREWORD 


CHAM of North America Incorporated has performed a Rocket Injector Anomaly 
Study under the NASA Contract NAS3-23352. ^ 

The purpose of the study was to modify, test and demonstrate a computer 
code for predicting three-dimensional two-phase spray flow and combustion 
in rocket engines. The modified computer code REFLAN3D-SPRAY (REactive 
FLow ANalyzer 3-D imensional, with two-phase spray) and results of parametri 
studies have been described in the following two volumes: 

Volume 1: Description of the Mathematical Model and 

Solution Procedure; and 
Volume 2: Results of Parametric Studies. 

Transfer of the code to NASA LeRC computer center, and preparation of a 
user's manual are recommended as next steps of the study. 

The authors wish to thank all those who have contributed to this work. 

In particular, thanks are due to Larry P. Cooper and Ken Davidian of the 
Communications and Propulsion Section of NASA LeRC; and to Laurence Keeton, 
Jack Keck, Kelli King, Janet Siersma, and Ronni Rossic of CHAM NA. 
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1 . 


SUMMARY 


The liquid fuelled rocket engine combustors consist of an injector plate 
and a thrust chamber. The injector plate consists of a number of propellant 
injectors which are designed to atomize the liquid jets of reactants 
and to promote intensive mixing between the vaporized components. 

Figure 1.1 shows the schematic of the rocket engine and injector plate, 
with L0X-RP1-L0X unlike triplet injectors, considered. The purpose of 
the study was to demonstrate an analytical capability to predict the 
effects of reactant injection non uniformities (injection anomalies) on 
heat transfer to the walls. For this purpose an existing three-dimensional 
single-phase flow and combustion computer code (REFLAN3D: - REactive FLow 
ANalyzer, 3-D imensional ) has been modified for simulating two-phase flows 
in liquid propellant rocket engines. The modified code is referred to as 
REFLAN3D-SPRAY. 

The numerical model in the final code assumes instantaneous evaporation of 
oxygen jets and treats fuel drops as discrete drops of given size spectrum. 

It accounts for the liquid fuel (kerosene) jet motion, evaporation, its 
interaction with the gaseous phase, and combustion. 

The coupling between the liquid jets and the gas-phase includes: 

- momentum interaction between the high speed liquid jets and 
gas phase ; 

- energy interaction between hot reacting gases and cold evaporating 
jets; and 

- mass transfer between the evaporating jets and gas phase. 

The velocity slip between the liquid drops and the reacting gases is of 
primary importance for accurate predictions of flow and heat-transfer 
characteristics. 

The Eul erian-Lagrangean approach for simulating spray flow, evaporation and 
combustion has been selected. A nonorthogonal body fitted coordinate 
system is employed for accurate simulation of combustor geometry and near- 
wall processes. 
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The physical models in the care are: the k ^ e model of turbulence, 

one and two step finite-rate reaction models, and six-flux radiation model. 

The numerical solution procedure employs a control -volume approach, with 
staggered grid and upwind-differencing practices. The solution scheme for 
Eulerian set of equations is a modified SIMPLE algorithm. These modifications 
include inproved equation solvers, under- relaxiation practices, and order 
of solving equations for different variables. For instance, velocities 
are solved by a point-by-point (Jacobi) method, while all other equations 
are solved simultaneously over the whole field. 

The code has been checked for both numerical and physical considerations. 
Results of test calculations and parametric studies as well as recommendations 
for further improvements and verifications are presented in Volume II. 
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2 . 


INTRODUCTION 


One of the significant ways in which the performance level of the rocket 
engines can be improved is by the use of optimal injector design, advanced 
materials and cooling concepts that allow a significant increase in thrust 
level. Figure 2.1 illustrates typical geometry of the rocket engine and 
selected injector element configurations. Injectors usually take the form 
of a perforated disk at the head of the rocket engine combustor and can 
vary in shape and dimension. Injection characteristics are improtant for 
proper mixing, evaporation, and efficient and stable combustion. The 
objective of the present project was to develop a predictive tool for rocket 
injection anomalies study. 

In the past, the design of injectors has been primarily based on experimental 
tests, experience and intuition. More recently, injector design capability 
is being improved further by employing computer codes capable of predicting 
two-phase flow, evaporation, chemical reaction and heat transfer within 
the complex geometries of liquid fueled rocket einges. Performance 
characteristics of liquid oxygen/hydrocarbons (LOX/HC) propelled engines 
can be studied by simulating combustion processes for different propellant 
combinations, injector elements, cooling systems, and pressure levels. 

A research program has been initiated by the National Aeronautics and 
Space Administration at the Lewis Research Center to provide a quantitatively 
accurate numerical modeling capability for the design and development of 
luqiid fueled rocket engines. The work has been performed at CHAM of 
North America Incorporated and resulted in a three-dimensional computer 
code (REFLAN3D-SPRAY) capable of predicting two-phase liquid fuel spray, 
combustion and heat transfer in engine combustors. An existing 3-D, 

CHAM computer code was modified to incorporate Eulerian-Lagrangean technique 
for two-phase spray simulation. 

Full details of the adopted mathematical formulation, physical models, boundary 
conditions and solution procedures are described in this report. An 
accompanying report, Volume 2, describes the results of various numerical 
consistency test cases, and parametric studies as well as the recommendation 
for further improvements and verification of the code. 
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Figure 2.1 Liquid Fueled Rocket Engine and Typical Injector Configurations 
















3. 


MATHEMATICAL BASIS 


3.1 Introducti on 

This chapter describes mathematical formulation of the two-phase flow, heat 
transfer and combustion phenomena handled by the REFLAN3D -SPRAY code. In the first 
section, basic principles of the Eulerian-Lagrangean (E-L) approach are presented. 
The second and third sections provide mathematical description of the physical 
processes in Eulerian and Lagrangean frames, respectively. The formulation is 
presented in fairly general terms so as to be applicable to the wide range of 
flow situations. 

3 . 2 Eulerian-Lagrangean Approach 

The mathematical formulation of the two-phase flow and combustion processes 
comprises the application of Eulerian conservation equations to the gas-phase 
and Lagrangean equations of droplet motion and thermal balance. The Eulerian 
part of the method involves solution of the following gas phase equations: 

- continuity equation; 

- three momentum equations; 

- energy equation; 

- turbulence kinetic energy and dissipation rate equations; 

- unburned fuel and CO mass fraction equations; 

- mixture fraction (composite fuel fraction) equation; and 

- radiation equations. 

The Lagrangean part is accomplished by integrating the droplet equations of: 

- motion; 

- heat transfer; and 

- mass balance 
along their trajectory. 

The spray combustion model assumes that the evaporating droplets act as 
spacially distributed sources of fuel vapor. The link between the phases 
involves mass, momentum and energy coupling and is mathematically expressed 
in terms of interphase transfer source terms in all Eulerian equations. 
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The complete solution algorithm for a gas-droplet spray flow is presented 
in Figure 3.1. 



Figure 3.1 Flow Chart of the E-L Solution Algorithm 

First grid and geometry data are specified and flow field variables are 
initialized. A droplet-free solution of the gas field is then obtained at 
the begining of the iteration cycle. This flow field is used for droplet 
trajectories, size and temperature calculations. The mass, momentum and 
energy interphase source terms are then determined. These source terms are 
now incorporated into the gas field equations. The new gas flow field solution 
is used again the the solution of the droplet equation. Thus there is a 
two-way coupling between the gas and liquid phase equations. The calculation 
porcess is continued until the converged solution is obtained. 
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3.3 Basic: Eulerian Equations 


3.3.1 Fluid Flow Equations 

The hydrodynamic equations are expressed in three-dimensional, cylindrical 
polar coordinates x-r-0, as follows: 


Continuity 

2p + JL (pu) +1 a(pvr) + 1 _i (pw) = jf, 

at 3x r 3r r 30 int 


(3.3-1)* 


x -Momentum 


_3 (pu) __3 (pu ) 1_ _3_ (rpvu) 1_ _3 (pwu) 

3t . 3x r 3r r 36 


3 ,p3u\ 1 3 / Hi3u \ 1 3 /Vi 3U\ 

3x l 3x 1 " r 3r 1 3r ' " r 30 V 30' 


= + JL + I JL + I _i 2 _1 / U /2L + I J.(rv) 1 3w x , } 

3x 3x ' 3x' r 3r ' 3x ' + r 30 ^ 3x' 3 3x m 3x r 3r * pK ' 


r 30 ' 


+ S 


int 


(3.3-2)** 


r- Momentum 


2L(pv) + 2L( puv ) +12 (rpv ) + 1_2 (p vw ) 
3t 3x r 3r r 30 


3 f V>3v \ 1 3_ / 3v \ 1 _ _2 fH. 3V\ 

3x v 3x ; " r 3r v M 3r' r 36 V 30' 


-3p , 3 /y3u\ 1 3 /ur3v\ 1 3 f /3w Wn 

3r 3x~ 11 3r' r 3r 1 3r ' r 36 m 3r " r' ; 


— ■ - (~ kk + ~) + pw 2 /r 
r r 30 r 


-2 3 r 1 3u 1 3 / \ , 1 3w\ , i-i , c v 

3 3r ^ 3x r 3r r 30 ^ ^ int 


(3.3-3) 


* See Nomenclature for explanation of symbols. 

**Note that for Cartesian coordinates, r-**>,3r = 3y and r3 = 3z. 
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0-momentum 


2 

3pw + _3 (puw) [ 1 9 (rpvw) 1_ _3_ (pw ) 3 / p3W x 1_ _9 , ry9w , I _i /P 9w\ 

3t 3x r 3r r 30 " 3x ' 3x' “ r 3r ' 3r ' ” r 90 V 30^ 


- liM + 

r 30 


3 /U 3u_\ 
9x V 30 ; 


1 3 


r 3r 1 V 30 


3v 


W 


D) + 


L_1 

r 36 


rH. 

V ^36 


+ 2v)} + H. (|S 


1_ 3v W' 

r 36 r 


pv w 

r 


1 9 rU / 3 U 

r 30 '3x" 


, 1 3 / » 

+ F3F < rv > 


L 3w 

r 30 


) + pk} + S 


w 

int 


(3.3-4) 


In the above equations, u, v and w are velocity components in the x, r and 
0 directions, respectively; p is the pressure; p and y are respectively the 
density and viscosity of the fluid mixture, which may be non-uniform. The 
normal-stress terms involving bulk viscosity (which is zero for most fluids 
including oxygen and RP1 fuel) have been omitted from the momentum equations. 


By multiplying equation (3.3-4) by r, and rearranging, a new 0 momentum 
equation can be obtained with angular velocity W = rw as a dependent 
variable. Experience indicates that some of the shear stresses expressed in 
the new form are easier to handle during numerical claculations. The 
transformed form of the G-momentum equation is: 

0-momentum 

M. + J. (puw) . I _i (prvW) 1 _3 (pw 2 ) J _, 3 Wx IJL /^.SMx 1 3 ,y 3 Wx 

3t 3x r 3r r 30 3x VP 9x' r 3r ' u 3r ; " r 30 V 30' 


- 1 M + 9 

“ " 7 30 + 37 





( u il) + 1 3 

' p 30 ' f 2 36 


/ 3Wx 2 3yv y 3v 2 3yW 

^30' R 30 r 36 " r 3r 


2 .1 _3 
3' r 30 



1 iiiv) 

r 3v 


VIt) + pk} + S 


w 

i nt 


(3.3-4a) 


Equation 3.3-4a is used in the code. For laminar flows, the velocity 
components are the instantaneous ones; and y is the molecular viscosity 
of the fluid mixture, For turbulent flows, it shall be assumed that the 
same equations are valid. For such situations, the time-mean values of all 
the flow variables and fluid properties are used, and y is 
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now the effective viscosity(which is the molecular viscosity augmented by the 
turbulent contribution). These comments also apply to the other equations 
which shall be considered in the following sections. 

Equations (3.3-1) to (3.3-4) govern the variation of the velocity components 
u, v, w, and the pressure p. In order to solve these equations, information 
about the variations of density p and viscosity y is required. This is 
discussed next. 

3.3.2 Density and Viscosity 

The density is related to pressure, temperature and the composition of the 
gas mixture through an equation of state: 

p = p (p, T, rcij ) (3.3-5) 

where T is the temperature of the fluid mixture, and the m.'s are the mass 
fractions of the component species of the mixture. Variations of T and m.'s 

vj 

are obtained either from part of the problem specification, or from the 
solution of additional differential equations. 

In REFLAN-3D, density is calculated from the following formula: 

P = |Hr (3.3-6) 

where ~ = Z Q (3.3-7) 

M j 

and where R is the gas constant and M- is molecular weight. Summation is 

>J 

taken over all species. 

For the evaluation of viscosity, laminar and turbulent flows have to be 
considered separately. For laminar flows, the viscosity is assumed to be 
a function of temperature and mixture composition: 

U = U (T, m.) (3.3-8) 

For turbulent flows, however, the problem is more complicated. The turbulent 
contribution (u t ) to the effective viscosity is a function of local quantities 
such as velocity gradients, etc. The evaluation of and the turbulence 
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model are discussed in the next section. 


3.4 The Turbulence Model 
3.4.1 The k[V€ Model 

The k 'v e model of turbulence is most commonly used for prediction of complex 
flow problems. The basis of the model is that it solves the transport equations 
for turbulence kinetic energy and turbulence dissipation rate. 

The turbulence model incorporated into the REFLAN3D code is the high Reynolds 
number k^e two-equation model recommended by Launder and Spalding (Reference 1)„ 
In the following section the governing equations are presented,, Details of 
the derivation can be found in the published literature (References 1, 2, 3) 
and are not provided in this report,, 

3.4 0 2 Governing Equations for the Ke Turbulence Model 

The basic differential equations for the turbulence kinetic energy k, and 
its dissipation rate e, are: 

JL (puk) + 1_ _9_ (rpvk) 1_ _9_ (pwk) _ _9 # r _ L _2 f r r ik\ 

9x r 9r r 99 “ 9x u k,eff 9x' r 9r v k,eff 9r' 

'FI ^ r k,eff F 9^ = G k “ pe (3.4-1) 

_2 (pue) 1_ _9_ (rpve) 1 3_ (pwe) _ , r . _ I _2 

9x r 9r r 99 9x ' ,1 £,eff 9x ’ r 3r e,eff 9r 

- i- ~ ^ r e,eff F 99^ ^ C l G k " C 2 pe ^ e/ k + peV ^ (3, .4-2) 

In the above equations is the generation term for the kinetic energy of 
turbulence and is given by: 
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- ro ,,3u% 2 . ,3V, 2 , ,1 3w , v 2 -, , ,3w , 1_ 3u, 2 ,3u_ 3v_, 2 

G k = y t {2{( 3^ } ^ 36 r } ^3x r 30^ ^3r 3x^ 


+ t 3w + I 3v _ Wx } 
v 3r r 39 r ; s 


(3.4-3) 


The quantities C D , Cj_ and C 2 are constants, r k gff and r^ eff are the effective* 
exchange coefficients for k and e, respectively, and are given by: 


r k,eff u eff /a k,eff 


(3.4-4) 


r e,eff y eff^ a e,eff 


(3.4-5) 


where and cr £ are the effective Schmidt numbers for k and e and are 

assumed to be constant . 


3.4.3 Turbulent and Effective Viscosities 

The turbulent viscosity y. is related to k and e, via: 

y t = C D pk2>/e (3.4-6) 

and, the effective viscosity is given by: 

y eff y t + y £ (3.4-7) 

where y^ is the laminar or molecular viscosity of the fluid. Often y^ is very 
large compared with y^, and y e ^ f can be taken equal to y t without introducing 
serious errors. The first practice (equation 3.4-7) is employed in the REFLAN3D code. 


Recommended values (Reference 3) for the constants appearing in the above 
equations are: 

C D = 0„09; 

C x = 1.43; 

C 2 =1.92; and 


subscri pt 


"eff" , to denote effective values, is used explicitly when 


coefficients that are related to turbulent flows exclusively are involved. 


otherwise non-subscripted symbols will be used. 
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0 k,eff = 0,9 °* 


(3.4-8) 


The Schmidt numbers a £ for the dissipation rate of turbulence is calculated 
from: 


(3.4-9) 


e ,eff 


(C 2 - Cj) c D 


where k is the Von Karman constant, A value of 0.4 is assumed for k. 


3 . 5 The Energy Balance Equation 
3 „ 5 . 1 The Stagnation Enthalpy Equation 
The stagnation enthalpy h, defined as: 

ft = C T + Em.H. + i- (u 2 + v 2 + w 2 ) 

P J 0 2 


is a dependent variable in the energy 

3 (putt) 1 3 (rpv?i) 1_ _3_ ( pwfV) 

3x r" 3r r 30 


1 3 /!> 1 3?u _ <•„ , $ 

F 30 ( h F 3F } - S ?( + S int 


transport equation: 

J3 3?k 1 3 ,rPjo 3?u 

3x n 3x " r 3r ' 3r ; 


(3.5-1) 


(3.5-2) 


In the above equations, C is the mixture specific heat, H. is the heat of 

r J 

combustion of the j-th species, S^ represents the sum of ^11 gas-phase source 
terms including thermal radiation (discussed below) and s!?^ represents 
interphase energy transfer source term. 


Equation (3.5-2) has been obtained with the assumption that the exchange 
coefficients for the transport of the mixture and for heat conduction are all 
equal at a point, although they may vary from point to point. 


It should be remarked that under certain circumstances the variation of the 

. 'jo 

specific enthalpy h can be obtained without solving the differential equation 
(3.5-2). For example, consider the incompressible flow of initially unmixed 
reactants in an enclosure with adiabatic, impermeable, non-catalytic walls 
and without any thermal radiation 0 A non-dimensional enthalpy can be 
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(3.5-3) 


defined as: 

^ o- 

h - h A 

^ = FTi" 

n F n A 


where the subscripts F and A refer to the inlet fuel and oxidant streams, 
respectively. The differential equation for the variable c|> h is then identical 
to that for the mixture fraction f in every respect, including the boundary 
conditions. Thus, ft and f are linearly related. If the reactants are pre- 
mixed and uniform in composition, then the enthalpy must be uniform throughout. 


Once the enthalpy and the species concentrations have been determined, the 
temperature T can be determined from equation (3.5-1) viz: 

ft - Em . H . - \ (u 2 + v 2 + w 2 ) 

T = (3.5-4) 


The specific heat C , however, is a function of temperature and the composition 

r 


c P - * m j y + v + y 2 * y 3 + y 4) 


(3.5-5) 


where aj, bj, cj, dj and ej are constants for each chemical specie. 

The Cp 'v T dependence is weak so the system (3.5-4) and (3.5-5) does not 
require iterations. Usually the previous iteratiorl value for T is used to 
calculate Cp and then (3.5-4) is used to calculate new temperature. 


3.6 Thermal Radiation 
3.6.1 Introduction 

There exist few numerical procedures for handling the radiation integro- 
differential equations. Of these the "zone method" of Hottel (Reference 4) 
and "Monte Carlo method" (Reference 5) are well known and tested. Recently 
introduced "beam tracing method" (Reference 6) at this stage of development 
is not suitable for engineering applications. All three procedures require 
significant quantities of computer time and/or storage. 


The iterative nature of reactive flow calculations requires simple and fast 
algorithms for radiation calculations. The most commonly used method for 
radiation simulation is the "flux model" originated by Schuster (Reference 7) 
and Hamaker (Reference 8) in astrophysics. Flux methods are based on the 
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use of some simplifying assumptions for the angular variation of the radiation 
intensity. This allows the exact integro-differential radiation transport 
equation to be reduced to a system of approximate differential equations in 
the space variable only. These are ideally suited to numerical solutions 
simultaneously with the flow equations. 

The six-flux method is described below for three-dimensional flow calculations. 


3.6.2 The Six-Flux Method 

For an absorbing-emitting grey medium in local thermodynamic equilibrium the 
radiation transport equation can be written as (Reference 9): 


ft • VI r = -al + af + 47 / Idft (3.6-1) 

4tt 

where I is the radiation intensity, is a unit vector representing direction 

4 

of radiating beam, a is absorption, s is a scattering coefficient, E = a.T , 
and a is the Bolzmari constant. Central to the assumption of a flux model is 
the assumed variation of the intensity with direction. Assumption of constant 
I in a quadrant centered along each (+) coordinate direction results in six 
differential equations for intensities (Reference 10). 


dF 


(rl) 


df(rd) 

dK 

dX 

dL 

dX 


I 

r d0 


I dN. 

r d0 


= r {-(a +s)I+~ + aE+|-(I + J + K + L+ M + N)} 
= r {(a + s ) I + — +aE-g-(I + J + K + L+ M + N)} 

= - (a + s)K + aE + g- ( I + J + K + L+ M + N) 

= (a + s)L - aE - | (I + J + K + L + M + N) 

= - (a + s)M + aE + |- (I + J + K + L+ M + N) 

= (a + s)N - aE - (I + J + K + L + M + N) 


(3.6-2) 


-} 

The I, J, K, L, M and N are the radiation fluxes in y, y , x , x ,0 , 0 
directions, respectively. 
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The net radiation heat fluxes, used in the energy equation source term can be 
evaluated from: 


V 


J = - 


a + s + 1/r 


d ( I + J) 
dr 


= l< - L 


2 d(K + L) 
a + s dx 


(3.6-3) 


Q 0 = M - N = - 


2 

a + s 


1 d(M + N) 
F dF — 


It can be seen that only composite fluxes (I + J), etc„ are used for the heat 
flux calculation. It is convenient, therefore, to add appropriate equations 
of the system (3.6-2) and arrive at the working equations for the " composite 
fluxes" R x , R and R z viz: 

H 1 dR v c 

( R x- E ) + !< 2 R x- R y- R z) ) 


_i (. 

dr ' 


a+s+p- 


1 _2 f 1 

r 39 ^a+s 


dR 




dr 


1 dR z, 


de^ 


r {a (R y - E) + f (2R y - R x - R z )} 
- {a (R z - E) + | (2R z - R x - R y ) } 


(3.6-4) 


where: 

R x h \ (I + J) 

Ry = \ (K + L) (3.6-5) 

R z = \ (M + N) 


The composite radiation-flux equations are easy to solve second order (diffusive 
type) ODE, They are coupled with the stagnation enthalpy equation through the 
presence of R x , R y and R z in the source term of the latter. Indirect coupling 
to other equations also occurs through the temperature appearing in the 
emissive power E and through any dependence of the coefficients a and s on the 
local quantities m^ u , m CQ , etc. 

The contribution of the radiation fluxes to the stagnation-enthalpy source 
term is given by: 


(Sft-) 


tr radiation 


= 2a {(R y - E) + (R x - E) + (R z - E)} 


(3.6-6) 
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Since information on the variations of the coefficients a and s with other 
quantities is often scarce and imprecise, they have simply been assumed to 
be constant in REFLAN3D. However, should knowledge of the variations be 
available, it is not a major problem to incorporate these into the code,. 

It may be noted that when the absorption coefficient a is very large, R , 

A 

Ry and R z become very nearly equal to the direct emission E. 

3.7 Chemistry Equations 
3.7.1 Reaction Models 

There are three levels of complexity for simulating the chemical kinetics in 
the REFLAN3D code: 

- instantaneous reaction assumption; 

- one step reaction mechanism; and 

- two step mechanism with intermediate CO. 

The rate of chemical reaction is assumed to be governed by an Arrhenius 
expression for laminar flows, and bly the "Eddy-Break-Up" model (Reference 11) 
for turbulent flows. 

Details of the combustion models and the chemistry equations are described 
below. 


3.7.2 Conservation Equation for a Chemical Species j 
The conservation equation for a chemical species j is: 


_3_ (pum.) 1 3_ (rpvm.) 1 3_ (pwm.) _3/ r Mil 

9x J r r J r 30 J “ 3x v j 3x J ' 


_ 1 — i. r rr M -j ) _ L 

r 3r ' j 9r 3 1 r 


9 

36 


(r . - 

J r 


1 3m 


36 


f j) = R. 


(3 o 7-1 ) 


where m. is the mass fraction of chemical species j; R. is the mass rate of 

J J 

creation of species j by chemical reaction, per unit volume; and r. is the 

vJ 

exchange coefficient. For a chemically-i nert species R. is, by definition, 

J 

zero. 
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3.7.3 Mixture Fraction and Inert Species Equations 


The simple single-step reaction of a pair of reactants (called here for 
convenience fuel (fu) and oxidant (ox)) can be expressed in the following 
stoichiometric relation: 

1 kg fu + s kg ox •+ (1 + s) kg products (3.7-2) 

where s represents the stoichiometric oxidant/fuel ratio and is a constant 
for a given reactant pair. It should be noted that the reaction equation 
(3.7-2) does not impose any restrictions on the constitution of the reactants; 
thus the reactants may be mixtures, e.g. fuel e CO + + Ng, oxidant = + Ng. 

An important consequence of the simple chemical reaction assumption is that 
the mass rates of creation by chemical reaction of fuel, oxidant, and product, 
R fu , R qx and Rp r , are related through: 

R fu - R ox' s - -V /(1 + S) <3 ' 7 - 3) 

These relations can be made use of to yield conservation equations that have 
zero source terms. 


It is further supposed that the exchange coefficients of fuel and oxidant are 
equal at each point in the flow, although they may still vary from point to 
point. It follows that equation (3.7-1) for oxygen (j'eox) can be divided by 
s and subtracted from the corresponding equation for fuel (jEfu). 


The result is: 


_JL (puy) . 1 JL (rpvy) 1_ _3_ (pwy) __3 (T3y) 
3x r 3r r 30 3x 3x 


i. _i (rr 8i) _ i_± 

r 9r u 9r ; r a6 


(r ~ = o 

u r 30' u * 


where y = m fu - m QX /s 


(3.7-4) 

(3.7-5) 


This is an equation having a single dependent variable, namely y, and no source 
term; the two reaction-rate terms have cancelled out. To make equation (3,7-4) 
more general, a normalized dependent variable f ("mixture fraction") defined as: 
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(3.7-6) 


f . o < y < 1 

Y P - Y A - 

can be introduced. The F and A indices represent conditions at fuel and air 
entry, respectively (Figure 3.7-1). 


(F) m fu,A' m ox,A 
m N2, A' 


(A) 


~j Inlets 

m fu, F' m ox,A — " 

m N2, A' 



Exit 


Figure 3.7-1 Specification of Inlet Composition at F- and A-inlets 


The mixture fraction transport equation is written as: 

9 (puf ) , 1 9 (rpvf ) 1 9_( pwf ) d_ (r 9ft 1 9 / r 9f\ 

9x r 9r r 90 " 9x ' f 9x' " r 8r ' j f 9r 


1 

r 


J 

90 


(r 


1 9f 
f r 90 


) = 0 


(3.7-7) 


The physical significance of the mixture fraction is that it represents the 
mass fraction of fuel in any form, i.e. fuel that has reacted and that has 
not. Thus, for chemically-i nert flows, the mixture fraction f and the mass 
fraction of unburnt fuel m- would be identical. 

Note that all atomic elements H, 0, C, N (inert species) are governed by 
the source free transport equation identical to that of mixture fraction 
(3.7-8). For the equal exchange coefficient, zero source terms and linearly 
related boundary values, all inert species concentrations at any location 
are linearly related to f. Stoichiometric relations provide the auxiliary 
1 inear relations. 
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3.7.4 Stoichiometric Relations 


The linear relation between the inert species <f> and the mixture fraction f can 
be expressed as: 


f - f a 

*F - *A " f F - f A 


(3.7-8) 


where <J> stands for m^, m Q , m^, m^,, etc., A and F indices represent A and F 
inlet, respectively. 


For simple one-step reaction there are five species participating in the mixture 
composition viz: (fu, 0 2 , N 2 , C0 2 , H 2 0). If two differential equations (f and 
m fu ) are solved, four additional algebraic equations can be obtained from 
relation (3.7-8). 


If the A-jet composition consists only of air (a 
F-jet of fuel (C n H m ) and say water vapor (m^ p, 
for atomic elements are: 


ox m ox + 
'V,0,F> 


( 1 -a ox ) V and 
the linear relations 


= 12 . n' f 

(3.7-9) 

^18 m H 2 0 F + m ' V fjT 

(3.7-10) 

1 6 

= a ox ■ < a ox ” 18 m H 2 0p) ~ 

(3.7-11) 


On the other hand, stoichiometric relations can be written as: 


m C ” 12 n ' m fu + 44 m C0 2 

m H " "■ m fu + ll m H 2 0 

m m 16 ,32 

m n = To Mu n + 7T7T 

0 ox 18 H 2 0 44 C0 2 

where n 1 = and m' s 


m 

12 n + m 


(3.7-12) 

(3.7-13) 

(3.7-14) 


Combining equation (3.7-9) and (3..7-12) m^ can be calculated. From equations 
(3.7-10) and (3.7-13) m„ n and from (3.7-10^ and (3.7-14) m can be determined. 

OX 
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(3.7-15) 


m H 2 0 p 

m H 2 0 = M9 m' + -pp— ) - 9 m' m fu 

m rn = 44 n'f - 44 n' m r (3.7-16) 

6u 2 fu 

a nx 

itIq = a - f (p* — + 32 n' + 8 in') + . (32 n 1 + 8 m 1 ) (3.7-17) 

More complex formulae can be derived for arbitrary A-inlet and F-inlet 
composition. 


Finally, since all the species mass fractions must add up to unity: 

m fu + “ox + m C0 2 + ™H 2 0 + ra N 2 = 1 

Thus m M can be obtained. 
n 2 


(3.7-18) 


The technique of relating mixture composition to the mixture fraction f can 
be expressed on two examples of combustion regimes: 

- diffusion controlled (instantaneous combustion); and 

- kinetically influenced reactions 0 


These are described in the following sections. 


3.7.5 Equations for Diffusion-Controlled Reactions 

Before any reaction can take place, fuel and oxidant must be brought into 
physical contact with each other. It is thus expected that the reactant 
process is affected by both the rate at which mixing of the two reactants 
takes place and the rate of the chemical reaction itself. In situations 
where the fuel and oxidant are not initially mixed, the mixing rate, being 
in general much slower than the reaction rate, even in turbulent flows, 
controls the reaction process. This leads to the assumption that 
thermodynami c equilibrium prevails everywhere. Thus, wherever fuel and 
oxidant are in contact, reaction will take place until one or the other has 
been completely consumed. Since chemical kinetics need not be considered 
under these circumstances, the complete chemical state can be determined by 
solving only one source-free conservation equation for the mixture fraction. 
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Once the f-equation (3.7-7) has been solved, the mass fractions of the other 
species can be obtained from algebraic equations. 


The mass fractions of fuel and oxygen are related to the mixture fraction 
according to: 

0 •= f - f st : m fu ■ 0 

V ■ (1 - f/t st> m ox, inlet (3 - 7 '- 19) 

1 > f i f st : m fu = - f/fst */' 1 - f st> 

m ox = 0 (3.7-20) 


In the above equations, f t is the stoichiometric value of the mixture 
fraction and is given by: 


m 


fu. 


f st = 1 7 { 1 + 5 <m > inlet } 

ox 


(3.7-21) 


The above equation is readily obtained from equation (3.7-6), noting that for 
stoichiometric mixtures, m fu = m ox /s, and consequently y ~ 0. The physical 
significance of f $t is that the locus of all points where f has the value 
f st defines the maximum reaction contour and hence the "flame envelope". 


The linear relationships between the mixture fraction and the various mass 
fractions can be used to calculate m^ , m^ g, etc. It should be remarked 
that these relationships are, strictly 2 spealing, valid only for instantaneous 
values. Thus, for turbulent flows, where only the time-average values of 
the mixture fractions are available, the relationships of Figure 3.7-2 can be 
used only for time-average mass fractions. This is assumed in REFLAN3D code. 


Some inaccuracy is involved in this assumption; it can be removed to some 
extent by, for example, solving an additional differential equation for the 
mean-square fluctuations of f (References 12, 13)„ 
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(A) Inlet ►. f (F) Inlet 


Figure 3 . 7-2 0 Diagram Showing the Variation of Mass Friction with f for a 
Diffusion-Controlled Reaction 

3.7.6 Equations for Kinetically-Infl uenced Reactions 

When the chemical process is not physically controlled (as, for example, in 

the case of premixed fuel and oxidant mixtures) the rate of the chemical 

reaction will be the influencing factor. Since the assumption of thermodynamic 

equilibrium can no longer be used, solution of a source-free equation alone 

is not sufficient to determine all the species mass fractions. In addition, 

a conservation equation for one of the species must be solved. This could be 

any one of the three: m f , m , or m 

J f u ox pr 

In REFLAN3D, for one-step reaction, two differential equations are solved: 
that for f, equation (3.7-7), and that for m^, equation (3.7-1) with j 5 fu. 
Under some circumstances, the solution of the f-equation may not be necessary 
even though the reaction is kinetically-controlled. For example, when a fuel- 
oxidant mixture of uniform composition is admitted into a chamber of impermeable, 
non-catalytic walls, under steady-state conditions, the source-free equation 
has the trivial solution: f is uniform everywhere and equal to the inlet value. 

Once the mixture fraction f and the fuel mass fraction m^ u have been determined, 

the mass fraction of 0 2 , C0 2 , H 2 0 and N 2 can be obtained from stoichiometric 
rel ations . 
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The variation of the species mass fractions with f for a kinetically-influenced 
reaction is shown in Figure 3.7-3. 



(A) Inlet ►- f (F) | n | e t 


Figure 3.7-3. Diagram Showing the Variation of Mass Fractions with f for 
a Kinetically Influenced Reaction 

NOTE: The dotted lines show the variations for a diffusion controlled 

reaction (as in Figure 3.7-2); the variation of m N is the same 
for both. 


3.7.7 Two-Step Reaction Model 

A general hydrocarbon ( C n H m ) fuel oxidation is employed in two-step parametric 
reaction scheme as follows: 

C n H m + + (1 - X) m j - 4 - - -) 0 2 ->■ An CO + (1 - X) n C0 2 + j H 2 0 

Xn CO + 0 2 -> Xn C0 2 (3.7-22 

where m and n represent C n H m -hydroc.arbon composition (e.g. for C 3 H g n=3, m=8), 
and A is the reaction scheme parameter (A=0 for single step chemical reaction, 
and A=1 for two step chemical reaction). 


For certain cases (e.g. for heavy liquid hydrocarbons) 0<A<1 can represent 
simultaneous homogeneous-heterogeneous reactions at the droplet boundary 
1 ayer. 
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There are six species (C^, C0 2 , CO, H.,0 and N 2 ) involved in the above reactions 
which have to be determined at each control cell of the combustion chamber. 


Three properties viz: mixture fraction f, defined as a "total fuel" (burned 

or unburned), mass fraction of unburned fuel i» c H , and mass fraction of 
carbon monoxide m c0 are obtained from the diffePeStial equations. Three 
additional equations are obtained from the stoichiometric relations: 

12 12 

m r = 12 n ' n V H + 2F m C0 + 44 m C0, 

0 n m 4 


m u “Hi nip M I 2 

H C n m + 18 m H 2 0 

16 m . 16 m , 32 
m 0 = m 0 2 + 28 m C0 + 18 m H 2 0 44 m C0 2 


(3.7-23) 


where n' = 


n 


and m 1 = 


m 


12 H~r¥ a,,u " TET+m 

3.7.7 Reaction Rate Expressions 

Two different options are considered for the reaction rate expression: 

- Arrhenius expression; and 

- Eddy-Break-Up model (EBU). 

The Arrhenius expression for the bi -molecular reactions can be written as: 


fu ,L 


= -pAM m , ( 


FU v M 


p m F!l a p m S 

FU ' ' 0, ‘.) exp (-E/RT) 


L ) (- 


(3.7-24) 


FU 


M 


OX 


where A is reaction rate constant, E is the activation energy, R is the 
gas constant, Mp^, are molecular weights and a and 6 are reaction order 
constants. Similar expression holds for m rn combustion rate. 


For turbulent flows, apart from the usual need to use time-averaged quantities 
and effective transport coefficients, the effect of turbulence on the reaction 
rate should be separately accounted for. The eddy-Break-Up model of 
Spalding (Reference 11) assumes that the gas is composed of alternating 
fragments of unburned fuel and almost fully burned gas (premixed flame). The 
chemical reaction is supposed to occur on the interfaces between these two 
inkds of gragments. 


The rate of reaction is supposed to depend upon the rate at which the fragments 
of unburnt gas are broken into still smaller fragments by the action of 
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turbulence, and is taken to be proportional to the rate of decay of turbulence 
energy. Thus, the rate of reaction is: 

R fu,EBU = - C R <=** (3.7-25) 

where C R is a constant, and g represents the local mean-square concentration 
fluctuations. 


The concentration fluctuation g can be determined by solving an additional 
equation having a form similar to that for the turbulence quantities k^e 
(References 12, 13). Magnussen (Reference 14) proposed an algebraic expression 
for g in the form: 


1/2 • r 

g = min (m 


m 


ox 


B m 


fu 1 


1 + s 


er) 


(3.7-26) 


where B = 4„5 is an empirical constant. 


For turbulent flows, the lower of the two reaction rates is often taken. 

3 . 8 General Form of the Governing Differential Equations 

All the differential equations discussed in the preceding sections are elliptic 
in nature and can be conveniently presented in the general form: 


sf (Wl) + F W (rW> + F55 


9 /T.ScK 
3F < 


l_ _ 
r 3r 


3 ,rl\ 3cf> 
~ ( ^ 3? 


1 


rW( r +fle) = s 




(3.8-1)* 


* I. As a reminder, note that for Cartesian coordinates: r-*», 3r^3y, and 

r30s3z„ 

2o The differential equations for the three composite radiation-fluxes 
are of one-dimensional form. Thus the r- and 0-derivatives of this 
equation are absent for R x , the x- and 0-derivatives are absent for 
R , and the x- and r-derivatives are absent for R . 
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Table 3.8-1. Summary of Equations for Three-Dimensional Flows 


EQUATION 
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In the above equation <f> identifies the dependent variable; B is identically 
equal to either the mixture density p or zero; is the appropriate exchange 
coefficient for the variable <j>; and is the source term which includes both 
the sources of <f> (positive or negative) and any other terms which cannot find 
a place on the left-hand side of the equation. Table 3.8-1 summarizes the 
equations in the form they are solved for,, Some notes about these equations 
now follow: 

- For laminar flows, instantaneous values of flow variables and 
molecular values of the exchange coefficients are used. 

- For turbulent flows, time-averaging values of flow variables 

and effective values of the exchange coefficients must be employed. 

The interphase transfer source terms will be discussed in chapter 4 0 6 0 

3 . 9 Basic Lagrangean Equations 
3.9.1 Introduction 

Drop life histories must be calculated in order to determine heat, mass, and 
momentum transfer. This is particularly important for pressure atomized 
injectors, where a significant proportion of the initial momentum in the flow 
is carried by the liquid phase and is transferred to the gas phase only by the 
drag force on drops. Since spray calculations are complex, the computation of 
droplet characteristics is represented in a relatively simple model. Droplets 
are assumed to be spherical and non-deformative ,with uniform conditions within 
each droplet volume. The droplets are divided into a number of drop size ranges 
and a system of differential equations is solved for each range. 

This chapter presents basic equations of droplet motion, heat and mass transfer. 
3.9 0 2 Droplet Distribution Model 

The present model assumes that the fuel is injected into the combustion 
chamber as a fully atomized spray which consists of spherical droplets. 

The droplet-size distribution within the spray is represented by a finite 
number of droplet parcels of specific droplet diameter. At the atomization 
point droplet initial sizes, velocities and temperatures are specified; these 
are subsequently tracked in a Lagrangean fashion as they traverse, heat-up 
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and vaporize within the combustion chamber. 


The number of droplets in each parcel, , is calculated according to a 
specified distribution function, that relates the probability Pj of finding 
a droplet of diameter d. to the droplet diameter: 

\J 

6m, P. 

N d = TET ( 3 - W) 

l P. p. IT <(3/6 

j=l J 3 

where NS is the number of droplet sizes and dnn is the mass of fuel introduced 
at the injection location. 


The most commonly used distribution function is a Rosin-Rammler function 
(Reference 15) (Figure 3.9-1) defined as: 

P R = a ^ (&) a 1 exp (- (?/) (3.9-2) 

K D 0 D 


where a is an 
diameter. 


Figure 3.9-1. 


empirical parameter (typically 1.5 < a < 3) and D is the main droplet 



The above consideration implies that the droplet trajectory, as well as heat 
and mass transfer within the spray, can be determined from the solution of a 
set of ordinary differential equations describing the behavior of each droplet 
parcel; these are given below. 
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3.9.3 Droplet Trajectory 

For each droplet size the momentum balance equations are written as (References 
16, 17): 

du d Ad 

m d Tt ' C D P T (u - u d> 1“ - “dl 


dv . /Lj Wj 

m d Tf' C 0 OT < v - v d> l v - v dl + "d T 


(3.9-3) 


dw d A d 

m d "dir = C D p T ( w " w a ) I w - I + — w r 


7 d 1 d r d 


Where Cp is the drag coefficient, u, v, w are gas velocity components, u^, v^, 
are droplet velocity components, is the droplet cross-section area , 
m^ is the droplet mass and r is the droplet raidal position. 

The other terms contributing to aerodynamic forces on the droplet include: 

- pressure gradient; 

- Magnuss force; 

- Saffman lift force; 

- Basset force; and 

- gravity forces, etc. 

These are all neglected because they are of the order of gas/droplet density ratio. 


The drag coefficients Cp depends primarily on the "relative" Reynolds number: 


Re 


P . D . U - U 


(3.9-4) 


where D is the droplet diameter. For evaporating droplets ,Cp can be calculated 
from the formula (References 18,19): 


C D = |t U - °- 15 Re°° 687 ) / (1 + 13) 


(3.9-5) 


for Re up to 1000. The Spalding (transfer) number B, given by: 

B=C y ^I (3.9-6) 
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can be significant for burning solid fuel particles. For evaporating droplets B 
is close to unity. In the above equation, C v is the specific heat of diffusing 
vapor, AT is temperature difference, and Q L is the latent heat of vaporization. 


If droplet velocities are established, the droplet trajectory (x^, y^, z^) can 
be obtained from simple relations: 


dx d dy d dr 6 

dt “ u d’ "dt “ V dt " w d 


(3.9-7) 


3.9.4 Droplet Heat Transfer Equation 

The vaporization process for a droplet moving in a high temperature gas 
stream is described in terms of two regimes: 

- heat-up period with raising droplet temperature, T^; and 

- equilibrium vaporization period with constant T^. 


A sketch of the droplet evaporation process is shown in Figure 3.9-2. 



Figure 3.9-2. Sketch of the Droplet Heat and Mass Transfer During the 
Vaporization Process 


At typical injection temperatures, the fuel concentration at the droplet 
surface is low and there is negligible mass transfer between the phases. 

As the liquid temperature rises the rate of mass transfer rises too, with the 
maximum liquid temperature occuring at the surface. Later in the process, due 
to the heat absorption for vaporization process, the droplet temperature reaches 
the so called "wet bulb temperature" which is almost uniform within the droplet 
mass. 
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At this stage intensive evaporation takes place with the temperature remaining 
constant at T^. These two regimes are accounted for in the mathematical model 
of the drop heat and mass transfer processes. 


Assuming quasi-steady state conditions, uniform conditions within the 
droplet volume and spherical shape of the droplet, the heat transfer process 
can be described by the following equation (References 20, 21, 22). 


i^d 

dt 


° 


(Nu . K (T - T d ) - aw (T^,, - T*) - Qj_) 


(3.9-8) 


where Nu = 2 + .06 Re 0,5 Pr 0,33 



d m d 
dt 


(3.9-9) 


(3.9-10) 


Here, T d is the droplet temperature; K, C d and L are thermal conductivity, 
specific heat of the liquid and latent heat of vaporization, respectively; 

represents the rate of the energy transfer with the evaporating mass, 
a, a and co are the Bolzman constant, droplet absorptivity and the view factor; 
T wall t * ie avera 9 e wa11 temperature. The 6/p d D factor represents droplet 
area to mass ratio. 

The rate of mass transfer dnip/dt is calculated from the droplet evaporation 
model which iis described in the next section. 

3.9.5 Droplet Evaporation Model 

Many single droplet and spray evaporation models exist ( see for instance 
References 23, 24, 25). The droplet dimunition rate is conventionally 
expressed by: 

d m , 

= tt D . K (3.9-11) 

where D is the droplet diameter and K is the evaporation constant (Reference 
25 ) o 

For computational purposes, equation (3.9-11) can be rewritten in terms of 
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rate of decrease of droplet diameter: 


d (D ) 

~dt = PD” K 


(3.9-12) 


The basic formula for the evaporation rate constant is expressed as (References 
25, 25): 

K = Sh . (p v K d ) (m y - m voo ) (3.9-13) 


where Sh is the Sherwood number, K Q is the diffusion coefficient, m y the mass 
fraction of vapor at the droplet surface, and m vro the mass fraction of vapor 
in the free stream. The Sherwood number is calculated from the Ranz and 
Marshal (Reference 27) formula: 

Sh = 2 + .6 Re 0,5 Sc 0,33 (3.9-14) 


For the high transfer rates Sh is modified by the convection factor and is 
expressed as: 

Sh = ;(2 + .6 Re° o5 Sc 0,33 ) .. (3.9-15) 


where B is the Spalding (transfer) number, given by: 

n c d( T -W 

B = c 


(3.9-16) 


and where C d is the vapor specific heat and T $at 
temperature. 


i s 


the wet bulb absolute 


3.9.6 Droplet-Wall Interaction 

The treatment of droplets impinging on the combustion chamber walls is one 
of the important difficulties in modeling spray flames. If the droplet hits 
the solid wall, a number of possibilities exist, e.g. the droplet may shatter 
into small ones which become re-entrained or it may adhere in the form of 
the sphere or a thin film, which subsequently evaporates. In the present 
method it is assumed that droplets adhere at the point of impact in spherical 
form, and that the heat and mass transfer processes continue to obey equations 
(3.9-8) and (3.9-12). It is recognized however that refinement of this approach 
will ultimately be necessary. 
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4. THE NUMERICAL SOLUTION PROCEDURE 

4 • 1 Introduction 

The system of coupled nonlinear partial differential equations (PDE) and 
ordinary differential equations (ODE) described in the preceding chapter can 
only be solved by numerical methods. 

This section presents numerical techniques employed for solution of both PDE 
(Eulerian equations) and ODE (Lagrangean equations) types. 

The finite-difference technique and a modified form of SIMPLE algorithm (Reference 
28, 29) are used for solution of the gas phase Eulerian equations. The resulting 
system of algebraic equations is then solved by efficient equation 
sol vers. 

The ODE equations describing the droplet behavior are first integrated 
analytically with the quasi -steady state assumption. The integration process 
is carried for each droplet class separately along its trajectory. During 
this process interphase transfer source terms for the gas-phase equations 
are calculated. This section describes both PDE and ODE solution practices 
in a detailed manner. First, however, geometry representation and space 
subdivision practice are presented. 

4 . 2 Geometry and Computational Grid 
4.2.1 Coordinates and Grid Lines 

The REFLAN3D code can handle both cartesian and cyl indrical -polar coordinates 
in both orthogonal and non-orthogonal form„ A non-orthogonal ity can be 
simulated only in the radial direction on the x-y (x-r) plane,, The code has 
also built-in moving grid options for simulating reciprocal processes 
(compressors, I.C. engines, etc.). 

The coordinate systems considered are: 

- x-y-z cartesian coordinates (Figure 4.1); and 

- x-r-8 cylindrical polar coordinates (Figure 4.2). 
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For cylindrical coordinates y = 0 does not necessarily coincide with the axis 
of symmetry, i.e. the radius r and radial distance y at any point may differ 
by a constant. 

The finite difference grid is formed by three sets of surfaces; perpendicular 
to coordinate directions. The nonorthogonal grid is formed by "piece-wise" 
surfaces combining r / r max points in the x-y (x-r) plane (Figure 4.3). The 
grid distribution can be nonuniform in each coordinate direction. 

4.2.2 The Staggered Grid Practice 

The REFLAN3D code employs the so-called staggered grid practice (References 
30, 31), in which all variables except velocities are calculated at the grid 
nodes (centers of the control volumes) (Figure 4.4). The velocity components 
are calculated at the cell faces of the control volume (designated by arrows 
in the diagram). A "backward boomerang" arrangement in the code implies that 
the velocities placed at west, south and low cell faces are assigned to each 
i, j, k node. The appropriate grid cell areas are also staggered in a 
"backward boomerang" fashion and for each node (i,j,k) calculated as A^, A<- and A^ 
for x, y and z, respectively (Figure 4.4). 

The derivation of the finite-difference equations will be illustrated in the 
first instance for a scalar variable, using the "finite volume" approach 
employed by Spalding (Reference 32). The porosity technique, used for 
representing complex geometries, allows any control volume to be fully open, 
partially open or. fully blocked. The porosity concept is presented in the 
following section. 

4.2.3 Porosity Concept 

In many engineering problems the boundaries of the domain are irregular. Also, 

there can be internal flow obstacles. In most of the known approaches 

complex curvilinear orthogonal body-fitted coordinates are used (References 

33, 34, 35) or nonorthogonal coordinates for simple geometry configurations 

are generated (Reference 36). The "finite domain" technique employs the so called 

"porosity concept" for simulating complex geometries. In this method every 

sub-domain is characterized by a set of volume (8) and cell face area (3 fl ) 

V H 
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fractions (cell porosities are usually: 0 £ 3 £ 1) that are available to the 

fluid flow. In Figure 4.5 a typical calculation domain of the flow over an 
obstacle is presented. In the present case finite porosities 0 £ 6 < 1 are 
specified within and around an obstacle. 

6=1 for fully open volume 

0 < 6 V < 1 for partially blocked volume 

6=0 for fully blocked volume 

A sample of a control volume within the blocked region is presented in Figure 4.5a. 
The volume 6 v V is used for gas flow calculations. A similar rule applies 
for the area porosities 6 AW , 6 AS , 6 al * 

4.3 Finite Difference Equations 
4.3.1 Motive of the Method 

The Finite Difference Equations (FDE) are obtained by integrating the Partial 
Differential Equations (PDE) over the finite volume (grid cell) and, for the 
transient equations, over the finite time interval. The following practices 
have been adopted for the spatial integration of the PDE's. 

- The flow variables stored at grid points are assumed to have 
stepwise profiles. 

- The "Upwind Differencing" (UD) (References 28, 29, 37) practice 
is employed for integrating convective terms. This implies that 
the scalar flow property required at the cell face is taken equal 
to that at the upstream grid point. 

- The cell-face velocity is considered as a cell-face average. 

No velocity interpolation is required for the calculation of each 
mass flux across the cell face. 

- Fully conservative formulation is used for the integration of each 
quasilinear PDE. This implies that all variables wlithin the 
iteration cycle are solved with the same (continuity preserving) 
fluxes. 
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For the time integration of the PDE's a fully implicit formulation is employed* 
i.e. the values of flow variables are taken to be those which prevail at the 
current time step. 


A seven-node finite-difference relation will now be derived connecting the 
value of a dependent variable <b at the node P (Figure 4.6) with those at the E, 
W, N, S, H and L neighbor nodes. Integration of the PDE will result in a 
linear formula of the form: 

a p^p = + a ^ ( ! ) i A | + a N^N + a S^S + a H^H + a L^L + (4.3-1) 


where a^, a^ s ..., etc. are called "link coefficients", SU and SP are the 
linearized sources and are given by: 


a P “ a E + a W + a N + a S + a H + a L“ SP c|) 


(4.3-2) 


is called the main (diagonal) coefficient. All the a's, SU and SP are treated 
as constants. 




Figure 4.6. Control Volume Notation for Cartesian and Cylindrical Polar 
Coordinates 
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There will be a set of equations like (4.3-1), each with individual coefficients, 
for every grid point and for each of the dependent variables. 


For a single variable a system of algebraic equations will be created which, 
in the matrix form, can be expressed as: 

t=S (4.3-3) 

The system can be illustrated on an example 5:3:4 (L:M:N) grid (Figure 4.7) , viz: 



Figure 4.7. Interpretation of the Storage Allocation for Three-Dimensional 
Computational Domain 

The corresponding matrix A to Figure 4.7 is presented in Figure 4.8. Each dot 
represents non-zero link coefficients. Diagonal lines filled with squares present 
link coefficients for the Z-cyclic (periodic) boundary conditions. 


Detailed discussion of the matrix inversion technique is discussed in 
Chapter 4.5, First, however, it is necessary to obtain expressions for the 
finite difference link coefficients (elements of matrix A in Equation 4.3-3) 
by integrating general transport equation in the form: 

+ af < pu<|)) + F3F (rpV((,) + r" a! (pwij>) 

' " -V ' 

convection terms 


<L/ r .3& 1 1 ^ 1 j fr 13t).o 

3x 9 3x' ~ r 3r 3r' “ r 39" 1 ^ r 39', “ <p 

y * 

diffusion terms 


(4.3-4) 
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source 

term 






Figure 4.8 Structure of the Influence Matrix for the Three-Dimensional Calculation 
Domain (see Figure 4.7) 

- Monzero Elements 

— □ — cj—tj — o — o — Cyclic Boundary Condition Links 
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The integration is discussed in the next section. 


4.3.2 Integration Over a Control Volume 

In order to present the derivations of the finite-difference equations s first, 
the x-direction related terms (equation 4.3-5) are discussed. 


<M. + _i 

3t 3x 


(pu^) 


— (p M.) + s 
3x K 9x ' 


(4.3-5) 


Then, a generalization to three-dimensional form is provided. Figure 4.9 shows 
a grid node P and its x-direction neighbors W and E. 



Figure 4.9. Control Volume and Notation for a One-Dimensional Transport 
Equation Integration 


The rectangular region represents the control volume used for the integration. 

At the w- and e-cell faces, convective fluxes C x and C x are defined as: 

w e 

C w = (puA) w and C e = (puA) e (4.3-6) 


where p is an "upwind" density: 


p w : 

if u 

.Pp* 

if u 


and the A's represent cell face area, i.e. 

AyAz for cartesian coordinates 

rArAG for cylindrical-polar coordinates 



(4.3-7) 


(4.3-8) 
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Integration of equation (4.3-5) in time and once in space yields: 


+ (pu4»A) e - (pu0A) w + 


- <‘ A - r 4, - < fl - r * H» - s 4, • mL 

T e Y w Y 


(4.3-9) 


Derivatives in the diffusion terms (in curly brackets) can be replaced by 
the appropriate finite-difference formulae: 


< AI 4, S7> “ A e r *e T 1 = °e <*E ' V 
e t 


(4.3-10) 


(AT |&) a A F (p - P w — = D (<f>p - <p. ) 
<b 3x w Y w 6X n w ' Y P Y W' 

w P 


(4.3-11) 


where D = A r/6X r and D w = A, T, /6X,, are called "diffusion link coefficients", 
e e <pe t w w cpw w 

The T^ and are calculated from appropriate linear interpolations 
between the node point values. The source term , whenever possible, is 
linearized in the form: 


s * • TOL ■ s % + sp 4, • * 


(4.3-12) 


The upwind differencing practice is employed for approximating convective 
terms : 


(pu<j>A) e - C e . 


(P^Aiw “ C w ‘ 


> P ifC e >0 


<jv if C < 0 
E e 

1f c w < 0 


cf),. if C > 0 
W w — 


(4.3-13) 


which can be rewritten using the following notation: 


fa,b] = max (a,b) 

as ( P u0A) e - fo, Cj <j) p + [0, -C e ] ^ 


and (pu4»A) ul - j0, C w | + |0, -C UI | <f > 


W JJ 


(4.3-14) 


(4.3-15) 
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Finally, equation (4.3-9) becomes: 

% - »V P + l0, Cj 4>p + do, -C e l 

- D e (* E - 4> p ) + D w (<j> p - <D W ) = SU 
pVOL 


where M = 
and M°= 


At 

p°VOL° 


At 


! <t> 
<t> 


E - f0. Cj 4>„ 

+ SP^ . 4> p 


fo, -C W U 4» p 


(4.3-16) 


(4.3-17) 


The integrated equation (4.3-16), analogous to equation 4.3-1, can be written 


as: 





Vp =: 

a^E 

+ a i/w + 

SU^ + M°4»g 

(4.3-18) 

where 

a P = 

M + a E + 

a W " SP (J) 

(4.3-19) 

and 

a E = 

D e + [0, 

-c e I 

(4.3-20) 


a W " 

D„ 1- f0. 

C J 

(4.3-21) 


Superscript "o" indicates "old time step" value. 


The above one-dimensional treatment of the differential equation (4.3-9) 
can now be generalized to three-dimensional flows. Consider the finite- 
difference equation: 

~ * a s^S + a H^H + a L^L + ^ (4,3-22) 

The various a.^ coefficients are given by: 

a P = a E + a w + a N + a s + a H + a L + M - SP^ 

a E " D e + - C el 

a W = D w + do, C w ] 

a N " D n + 

a s = D s + fib, C $ ] 

a H = D h + l 0, “ C h^ 

a L = ^1 + t® ; » "^l 


(4.3-23) 


(4.3-24) 
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The convective C's and diffusive D’s fluxes for the y and z directions are 
calculated from the expressions similar to that of equations (4.3-11) and (4.3-6). 


4.3.3 Modifications for Non-Orthogonality of Grid 

Figure 4.10a shows an example of the computational domain with nonorthogonal 
grid at the downstream part of the chamber. The non-orthogonality is allowed 
only in the axial (x-r) plane. Figure 4.10b depicts a selected control volume 
of the grid in the non-orthogonal region. 


The inclination of the cell face from the horizontal direction is denoted by 
an angle a. a varies with radius such that a = 0 at the axis of the chamber, 
and a = ex „ at the north wall of the chamber. 

The net convective flux of fluid across the inclined south face of the control 
cell shown in Figure 4„10 is: 


C s = p s A s 


p s A sx 


(4.3-25) 


where A g = horizontal projection of the area of the inclined north face; and 
A = vertical projection of the inclined north face area. 

SX 


A s and A gx are related to each other such that: 

A sx = A s * tan a (4.3-26) 

The axial velocity at the south cell face u g is calculated from a linear 
interpolation between the u p , u^, u s and u sp velocities,, Similar expressions 
can be derived for the north face. Note that in the orthogonal regions, A 

S X 

is zero and therefore equation (4.3-25) reduces to the standard convective 
flux. 

c s ■ \ v p < 4 - 3 - 27 > 
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4.3,4 The Momentum Equations 

The momentum equations have a finite-difference form similar to that of the 
general <j>-equation (4,3-20); but two main differences exist. First, the 
velocity components are at "staggered" grid locations ; therefore, the control 
volumes used for the velocity components are different from those for the 
other dependent variables. Second, the pressure-gradient term which forms a 
source of momentum is given a special treatment. 

Figure 4.11 presents a typical control volume arrangement for the u-velocity 
equation integration. 



NW , 

c y 

i NW 

— N, 

c Y 
i C N 



x i 

c w | 

. w 

! 

_ J 

i 

i p J 

c x 

.E 


• 1 
i 

o 

l£< 

1 

1 

1 

' ^>0.1 
a j 

I 

! 
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» 

s, 

\ 

U - Con 

itrot Volume 


Figure 4.11. Control Volume Arrangement for Integration of u-Momentum 

Note that the pressure node is always in the center of the <}>-control volume. 
However, the velocity for a nonuniform grid is not necessarily placed in the 
middle of the u-control volume. 


The finite difference equation form of the x-momentum equation is, with the 
pressure gradient term written out separately from other momentum sources: 


a P U P " a E U E + a W U W + a N U N + a S U S * a H U H + a L U L + SU u " A w ^ P P P lP + ^ U P 


(4.3-28) 


If velocities are solved by Jacobi iteration, it is desirable to rearrange 
the above equation to the following form: 


a c u c + a..u., + a..u., + a c u c + a u u u + a, u. + SU + M^u° 

Up . _tl li-JU 4 —^ Dyp p . p u ) 

where : 


DU 


P~ a 


(4.3-29) 

(4.3-30) 


P,u 
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The link coefficients can be written as: 
ap ' a E + a y + a N + a s + a H + a L + I" - SP„ 

a E = .5 {(D e + D p ) + [ 0 ., -C* - C*1 } 

a „ ■ -5 «D p + D w ) + to, C x p + C*J} 

a N = ' 6 t(D n +D n«> + lo., - I + i° - ■ - c n u i > (4.3-31) 

a s - -5 «D S + D sw ) +10., C*] + do., C* 1 ) 
a R = .5 ( ( D h + D p ) + lo., - Cjl + l0., -C^l } 

a L = ,5 {(D ] + D p ) + [ 0 ., C* J + LO., -C* w 1 } 

The r and 0 momentum equations have similar finite difference form. 


Note that the convective parts of a^, a<~, and a^ link coefficients are 
calculated for both convective subfluxes separately (see Figure 4.11). In the 
alternative approach, first a net convective flux is calculated as a sum of 
two subfluxes and then upwind principle is employed, e.g„: 


= - 5 + »„«) + -C 1 } 


This practice is used in several existing and widely used codes , e.g. C0RA3, 
TEACH., C0M3D and STARPIC. 

Figure 4.12 presents comparison between the practices for a selected flow 
configuration (C^ + Cp < 0 and > 0) 



(a) REFLAN (LINKS WITH U N , Ug and U w ) (b) OTHER (LINKS WITH U N and U w ONLY) 

Figure 4.12. Convective Link Coefficients for U p - Momentum Equations 
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In recirculation regions, the practice employed in REFLAN3D retains convective 
links with both North and South neighbor velocities u^, u<.. Other practices 
maintain the link only with u^. 

The practice used in REFLAN3D code is based on a full conservation principle, 
and, by employing more links, it is also more accurate than the other practices 
mentioned earlier. A similar treatment yields the link coefficients for 
radial v and circumferential w velocity components. 

The momentum equations are solved by using a modified version of the SIMPLE 
algorithm (Reference 28). In SIMPLE, the momentum equations are first 
solved with a guessed pressure distribution, denoted by p*, to give a first 
approximation to the velocity fields, u*, v* and w* (the starred-velocity 
fields). These velocities are approximate because they do not in general 
satisfy the continuity equation. Corrections to the pressure and velocity 
fields are then obtained from the solution of pressure-correction equations 
which are derived from the continuity equations (see Section 4.3.5 below). 

These corrections are such that the resulting velocities reduce continuity 
errors. The iterative is continued until convergence. 


4.3.5 The Continuity Equation: Pressure-Correction Equation 

It has been mentioned above that the velocities obtained from the momentum 
equations do not satisfy the continuity equation, and therefore require 
correction. The correction of velocities and pressure is discussed in this 
section. 

The continuity equation is: 

3 “ (pu) + ~~ (pvr ) + (pw) = 0 (4.3-33) 
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It can be written in the finite difference form as: 


PpVOL 

At 


P d °V0L° 

_ -JbL— 

At 


+ C. 


C X + C y - C y + C* - C* = 0. 
w n s h ! 


(4.3-34) 


where (to remind) C g = (puA) e , ... 


(4,3-35) 


Again, to start with, a one-dimensional problem is considered. From the 
momentum equations, the following relationships are obtained: 


u e = u e * + DU e ( p E " p P (4,3-36)* 

and 

U w = U w + DU w (p P " P W } (4.3-37)* 

Here p' is the pressure-correction , and DU's are the pressure-difference 
coefficients which were calculated from equation (4,3-30) during the solution 
of the momentum equations. 


For compressible flows additional density correction is performed. The 
correction practice employed is based on the pressure density relation: 

P = P* + (|£) . P' (4.3-38) 

where (|£) is calculated from equation of state. 


Substitution of p (4.3-38) and u (4.3-36 and 37) into the mass conservation 
equation (4.3-34) results in: 


VOL 

At 


(p* 


+ p 1 ) 
9p P ; 


p 0 VOL 0 

IE - 

At 


+ (u e * + DU e (p £ - pf,))(p e + If P^) A e + 


- <V + (p‘ p - p’)) (p w * f£ pj,) A w + ... (4.3-39) 

** Equations~{4".3-36) and (4.3-37) are approximate forms of the momentum equations. 
The exact form would be: 

u e = u e* + DU e (p P “ P E } + l . a i (u i “ u i* } 

where i indicates summation over the neighboring nodes; and a. are the 
finite-difference link coefficients from equation (4.3-29). 
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If an upwind differencing practice is used for cell-face density calculations 
and ignoring terms of order the pressure correction equation can be written 

as follows: 


a p Pp =: a E u E + V W + a N u N + 
where a p = + a y + 

a W = P w A w DU w + I°*» 


a E = p e A e DU e + i0 ” 


a N = p n A n DV n + l °*’ 


a S " p s A s DV s + & 0 *’ 


a H ~ p h A h DW h + i0 -’ 


a L = P 1 A 1 DW 1 + “ °* 5 


a S U S + a H U H + a L U L + SU + M* - M° 

a E + a N + a S + a H + a L 
UjA w • OP/3P) w 

ujA e . Op/ap) e 

V n JA n . (3p/3p) n 

V S 1 A $ . (3p/3p) s 
W h lA h . (3p/3p) h 
W 1 ]A 1 . (3p/3p) 1 


(4.3-40) 


(4.3-41) 


The source term SU represents net mass flux imbalance of the u*, v* , 
field. Thus: 


SU 


w 


+ C; 




+ c 


h 


'1 


w*-velocity 

(4.3-42) 


The purpose of the pressure-correction equation is to reduce this mass source 
to zero. When these sources are everywhere zero, the solution is just p 1 = 0. 

After solving equation (4.3-40) the corrected velocities are obtained from 
equations (4.3-36) and (4.3-37). (Similar equations bold for the v- and w- 
velocities.) Density is corrected via equation (4.3-30), and pressure, from: 

p = p* + p 1 (4.3-43) 
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In the derivation of equation (4.3-40), the approximate form of the momentum 
equation has been used. However, no error is thereby introduced. In a 
converged solution for a steady-state problem, the mass-sources and the pressure 
corrections will have become zero (in reality, a relatively small number) in 
the final iteration so that the solution is independent of what actually went 
into the derivation of the pressure-correction equation. 

Note that after solving p 1 , correcting velocities and calculating new, 

r\j 

continuity obeying, convective fluxes the remaining equations, k, e, n ... , etc. 
can be solved in a fully conservative formulation. 


4.4 The Radiation-Flux Equations 


The differential equations for the three composite radiation fluxes R , R 

x y 

and R z , equation (3.6-4) is much simpler than the general equation (3.8-1). 
The equation for the x-direction flux R can be written as: 

A 


d 

dx 


dR 

< r '-df 


■) 


+ S = 0 


(4.4-1) 


where r = l/(a+s) 


(4.4-2) 


and S = a(R x - E) + s(R x _ R y )/2 


(4.4-3) 


The finite difference form of thds equation is obtained by integration over 

the x-direction length of the control volume. The integral of the source 

term S is expressed as (SU + SP . R ) as before, and the complete finite- 

xp 

difference equation can be written as: 


a P R xp = a W R x + a W R x + SU 
K e w 


where a p = a w + a £ - SP 


a., = r /ax„ 
W e e 


a,, = r /ax 
W w' w 


(4.4-4) 


(4.4-5) 


The distances 5x and 6x are defined in Figure 4.9. 


Similar equations can be obtained for the y-direction and z-direction fluxes 
R^ and R z , respectively. 
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4 . 5 Solution Method of the Complete Equation Set 

For the NX, NY, NZ grid and NVAR dependent variables (u, v, w, p, k, e, ft ...) 
to be solved for, the finite difference equation set can be written in the 
matrix form as: 

As = B (4.5-1) 

where s is the solution vector containing all the unknown velocities, pressures, 
enthalpies, concentrations , etc. , B is the "source" vector and A is the 
coefficient matrix of the form: 

Figure 4.8. Matrix Arrangement 
for the Complete 
Equation Set 


where the form of D^, D v , submatrices is as that on Figure 4.8. 

It is seen that the matrix A has a large and extremely sparse structure. For 
a typical 3-dimensional problem with 20 x 20 x 20 grid and 10 dependent 
variables (80 000 algebraic equations) matrix A has 64 million elements. Of 
these, less than 48,000 are non-zero. 

Utilization of any direct solvers, such as Gaussian elimination (Reference 40) 
or band solvers (Reference 41), even for much smaller systems is inpractical 
for the following reasons: 

a) prohibitive storage; 

b) such a solver would have to employ partial pivoting because of 
zeros along the diagonal; 

c) solution time of direct solvers is proportional to the number 
of unknowns to the power 2 to 3; and 

d) the equation system is nonlinear so that several matrix inversions 
would be required. 
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The only alternative solution algorithm is to solve the (4.3-3) equations 
separately for each dependent variable. Even in this case, solution of 
equation (4„3-3) for three-dimensional flows is a formidable task. 

It has been mentioned before that the momentum equations are solved by the 
Jacobi, point-by-point algorithm (References 40, 42) „ The remaining equations 
including p'-equation, are solved simultaneously by the so-called "Whole Field" 
solution process at all control cells within the calculation domain. 

The simultaneous solution practice employed in the REFLAN3D code is described 
in the following section. 


4.5.1 The Whole Field Solution Practice 

In last five years, the commonly used "Alternating Direction Implicit", ADI, 
method (References 38, 24) for the solution of a system of algebraic equations 
has been replaced by Whole Field Solvers (WFS). 

Spalding and his coworkers (Reference 60, 61) have devised and used a WFS 
which is similar to Stone's Strongly Implicit Procedure, SIP, (Reference 44) 
but completely free from any adjustable (convergence-promotor) parameter called 
a-parameter. Subsequently Gosman and others (Reference 45 to 48) have 
investigated effectiveness of the WFS. All of these solvers are itenative 
field traverse methods involving a forward march for assembling the solver 
coefficients and backward march for back substitution. However due to the 
inconvenience of indexing practices, the marching directions are not 
reversed or alternated, as a result these solvers are not symmetric. 

A new fully symmetric solver which is totally independent of the a-parameter 
has been recently developed at CHAM by Przekwas (Reference 49) and employed 
in the REFLAN3D code. 
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Figure 4.13 illustrates the grid node arrangement and the nomenclature for 
the WfS discussion. 


N 



Figure 4.13. 



Grid Node Nomenclature for the WFS Discussion 


With the i» j, k indices omitted the FD equation can be expressed as: 


V = a E^i +1 + a w^i-l + Vj+1 + a S (f, j-l + a H <J) k+l + a lA-l + SU 


(4.5-2) 


The solution algorithm of the z-symmetric 3-D WF solver for the system of the 
FD equations (4.5-2) can be summarized as follows. 


Forward March : i) calculate <j> independent modified coefficients (in increasing 
i , j , k order) 

E = a E /D 
N s a N /D 
H = a H /D 
L = a L /D 

D E a n - a., E . i - a c N . . 

P W i-l S j -1 

ii) calculate the modified right-hand-side (in the order as above). 



Backsubstitution : iii) solve. tridiagonal equations (in decreasing i j order). 

" L *k-1 + ^ " H0 k+1 = BB + E(f) i+1 + N(f, j+1 

using TDMA or CTDMA (see next section) along each k-line„ 

iv) return to ii until number of sweeps or convergence criteria has been reached. 

This method has an additional advantage of implicit treatment of the periodic 
boundary conditions in the z-direction. It is also generally faster than any 
of the earlier mentioned practices. 


In the following section a TDMA (Tri -Diagonal Matrix Algorithm) and CTDMA 
(Cyclic-TDMA) are discussed. The algebraic equation (4.5-5) along any z^line 
(k=l . NZ) can be expressed in general form as: 


Vk-l + Vk + c k*k+l “ D k 

where A k = L k 

B k = 1 

C k = H k 

D k = 8B k + Vi+l + Vj+l 

and can be used as a basis for the following discussion,. 


(4.5-6) 


(4.5-7) 


4.5.2 The TDMA Algorithm 

The solution algorithm consists of two steps: 

a) forward march: where the TDMA coefficients are calculated as follows: 

-C,, 

v 

5l k a k-l 

- ~ (4.5-8) 


a k ~ B. + A, a 


6 = °k " A k B k-l 
k B k + A k a |< _ 1 


k = 1, 2, ... NZ 
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b) back-substitution: where values are calculated as follows: 

^NZ = 3 NZ 

k = NZ-1 , NZ-2 .... 2, 1 (4.5-9) 

\ = a k *k+l + 6 k 
4.5.3 The CTDMA Algorithm 

The Cycl ic-TDMA algorithm is used to solve the system of equations with 
periodic (cyclic) boundary conditions. Figure 4,14 presents a grid arrangement 
with cyclic boundary conditions. 


N— 1 N 1 2 3 k-1 k k+1 N-1 N 1 2 



+ B i <j> y +C-, 0 2 ~ D i 






A N ^ N-1 + b N 0 N + C N 0 1 = D |\ 


A k 0 k-1 + B k 0 k + C k 0k+1 " D k 



Figure 4.14. Grid Arrangement and Equation Forms for Cyclic Boundary Conditions 

Note that the (k-1) neighbor of <p 1 is 4> NZ . Detailed description of cyclic 
boundary conditions is discussed in chapter 5. Here only the solution algorithm 
is outlined. Similarly as for TDM A the solution algorithm consists of two steps 

a) forward march - where the CTDMA coefficients are calculated as 
follows: 
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for k = 1 


aj = D^/B^ 


f = d nz 


h = C i /B i 

and 

G = C NZ 

(4.5-10) 

°i = V B i 


H = B nz 


for k = 2, .... NZ-1 




a k = (°k “ A k k-l^ /W x 


F k = F k-1 " G k-l a k-l 


S k " C k /W w 

and 

G k = " 3 k-l G k-1 

(4.5-11) 

a k = A k 0 a k-l /W k 


H k = H k-1 _ G k~l a k-l 



where a = - A k 3' k _ 1 (4.5-12) 

b) back-substitution - where <p values are calculated: 

, _ F NZ-1 ~ ^ G NZ-1 ~ A N Z^ • a NZ-l f/L c: 10) 

NZ H NZ-1 " (G NZ-1 + A NZ* (3 NZ-1 + a NZ-l ^ 


^k = a k " B k ^k+l = a k ^NZ k = NZ_1 » NZ ~ 2 ’ •••> 2 » 1 (4.5-14) 

4.5.4 Under-Rel axation 

The calculation procedure described above involves the solution of nonlinear 
differential equations expressed in the form of the linearized finite-difference 
equations. Therefore, an iterative procedure has to be applied to continuously 

update the coefficients until a converged solution is obtained. If the changes 
in the values of the variables from one iteration to the next are large, there 
is a possibility that convergence may not be achieved at alio To keep these 
changes sufficiently small, the dependent variables are suitably under-relaxed. 
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There are two under-relaxation practices employed: 

1. inertial under-relaxation of dependent variables (u, v, k, e ...); 
and 

2. direct under-relaxation for secondary variables (y, P). 


Inertial under-relaxation of the dependent variables is achieved by adding 
an inertial term viz: - cj>*) to the finite equation in the following 

manner: 


Za njn + SU+JLJ* 


Za„ 


SP + I 


<P 


(4.5-15) 


where suffix n denotes all cell neighbors, and supscript * denotes previous 
iteration value of <p. The "inertia term" I is calculated as: 


I 


- P VO L 



(4.5-16) 


where p is the fluid density, VOL is the grid cell volume and Atp^ is the 
"false" time step specified for each dependent variable <p. The main features 
of equation ( 4 0 5-1 5 ) can be summarized as: 

1. In a converged solution cj) = cj>*, and therefore the final solution of 
the finite difference equation is not affected by the magnitude of I. 

2. has the dimension of flow rate,, 
f : 4> 

3„ The smaller the value of At c , , the heavier is the under-relaxation. 

Fcp 

4. I<j>* and I are included in SU and SP components of the source term 
so that the general form of the equation remains as Equation 4.5-2. 


The direct urider-rel axation practice is implemented by calculating the under- 
relaxed (|) U -value as a weighted average of just-calculated <j>-value and previous 
iteration cj>*--value in the form: 

(f> U = cyfr + (1 - a<j) ) 4>* (4.5-17) 


where a ^ is the under-relaxation factor for <p variable and has the value 
between zero and one (a = 1 implies no under-relaxation). 
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4.5.5 Calculation of Residual Errors and Convergence Criteria 


The residual errors are calculated for the equations of all dependent variables 
at each control cell, in the following manner: 

£(j) p = Za n 4> n + SU - SP <|> p * - (Za n ) (j) p * (4.5-18) 


where e^p stands for the residual error in the equation of variable <j> at point 

P, and summation is taken over all link coefficients: N, E, W, S, H 

The above equation is obtained from equation 4.5-16 by transferring all terms 
to the right-hand-side and equating them to the residual error. Three major 
quantities can provide the information on the convergence of the solution 
obtained and can be printed by the REFLAN3D code at any iteration. These are: 


1„ Maximum residual error 

RESMAX = max (e ) (4.5-19) 

ijk Mjk 

2 0 Global residual error 

RESSUM = Z e. (4.5-20) 

ijk Mjk 

3. Global absolute residual error 

RESSUMABS = Z |e. | (4.5-21) 

ijk M'jk 


Additionally an absolute difference between two consecutive iteration 
cj)-values is calculated as: 


DIFMAX - max (* 1jk - * 1jk *> 


(4.5-22) 


which is searched for over the entire calculation domain. 


The solution is regarded as converged when all above quantities take values 
below their prescribed limits. Usually it requires reduction of two, three 
and sometimes four orders of magnitude of the residual before this occurs. 
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4 . 6 Integration of the Lagrangean Equations 
4.6.1 Integration Practice 

The Lagrangean equations for the droplet motion, heat and mass transfer (see 
Chapter 3.9) can be expressed in a general form of ordinary differential 
equation (ODE) as: 


d <j> 

dt 


f - '* , d ) 


+ F 


(4.6-11 


where <j> and <f> . represent gas-phase and droplet property, respectively. 

9 ^ 

The E and F coefficients in general Lagrangean equation (4.6-1) can be 
obtained from equations governing droplet momentum (3.9-3), energy (3.9-8) 
and mass (3.9-11) conservation. The solution practice for the Lagrangean 
and Eulerian ones is totally different. The PDE's of the Eulerian part 
are solved as a boundary value problem while the ODE's of the Lagrangean 
part represent an initial value problem. 

Several techniques for solving a system of coupled ODE's of the droplet behavior 
have been used, including integration formulae (Reference 50), second order 
predictor-corrector schemes (References 51, 52), and fourth order Runge-Kutta (R-K) 
scheme (Reference 53, 54) and analitical integration techniques (Reference 20, 21, 22) 


In this study equations (4„5-l) are ;solved analytically. This technique offers 
considerable economy in computing time in comparison with the iterative fourth 
order R-K scheme (Reference 53). 


Integration of equation (4.6-1), assuming the gas property <j> is constant over 
the time of integration, yields: 

<f> = 4>g - (<f>g - <J> d ) exp (- — ) + E F (1 - exp (- —-) ) (4.6-2) 

where At is the time interval. It is worthy to note that the E factor represents 

a characteristic time scale (or relaxation time) for the process. If, during 

the integration process, the time scale ratio At/E is larger than say 20, the 

-9 

exp (-At/E) is smaller than 2.10 and the integration effort of equation 
(4.6-2) is greatly simplified. 


4-28 



4.6.2 Droplet Trajectory 

After determining droplet velocity at time t + At, the droplet position at 
time t + At is determined from: 


- ■■ o / i”i .1 o» At 

x d = x d ( U d ‘ U d ^ 2 


(4,6-3) 


where x d ° is the droplet position at the beginning of the time increment and 
U^ 0 is its initial velocity (at x^ 0 ). 


Figure 4.6-1 depicts consecutive droplet positions along its trajectory crossing 
the orthogonal Eulerian grid cells „ The integration process starts at the cell 
surface (position x^ in Figure 4.6-1) where cell droplet inlet properties 


U, , V . , W . , T, , D, , are solved. 

U IN a IN a IN IN a IN 


• N , NE 



Figure 4.6-1. Droplet Trajectory Passing the Computational Cell 


Next the estimate for the time interval for the integration process 
calculated as: 


,<U p + (V p + V„)/2 U d V dA 

At a INT o min ( Ax , Ay » Ax ’ Ay' 


i s 


(4.6-4) 


where a^y < 1 is the parameter approximately specifying the number of 
integration interval before the droplet reaches the exit position at the grid 
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cell face. As an example, in uniform gas and droplet velocity field, = 1/5 
would imply five time steps within the grid cell. 

Then the new droplet velocities, temperature and diameter are calculated from 
equation 4.6-1 and new droplet position is established from equation 4.6-3, 
(point 2 in Figure 4,6-11). Next, integration time interval At is calculated 
again and the integration process is repeated until droplet crosses the cell 
face (point 4 in Figure 4.6-1). 

At this stage linear interpolation is employed between points 3 and 4 to find 
droplet properties (u^, v^ , w^ , T^, D^) at the cell face (point 4 in Figure 
4.6-1). A droplet exit point "Ex" has been found and the integration procedure 
is completed by calculating interphase sources between the liquid and gas phase. 
The mass transfer, for example, is calculated as: 

iTi INT = 6 p d (°dIN ' D dEX ^ rt d (4.6-5) 

where is the number of droplets traversing the grid cell in unit time. 

Detailed discussion of the interphase source calculation is discussed in the 
following section. 

The exit droplet position in the P grid cell (point 4 in Figure 4.6-1) is 
considered as an entry point for E-grid cell and the integration process 
starts again. The droplet trajectory is traced until the first of the following 
conditions is met, 

a) droplet diameter diminishes to zero; 

b) droplet leaves through the exit of the calculation domain; or 

c) droplet hits the chamber wall where it evaporates. 

A special treatment is required for tracking the droplet on a nonorthogonal 
grid. Figure 4.6-2 depicts a typical droplet trajectory and its inlet (IN) 
and exit (EX) cell boundary intersection locations,, 
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Figure 4.6-2. Droplet Tracking on Nonorthogonal Grid 

In the case of orthogonal coordinates the cell boundaries XM p , XM p , YM^, YM p 
(Figure 4.6-2) were constant within the grid cell and independent on the 
droplet location. For nonorthogonal coordinates, however, at every integration 
step a continuous update of YM p and YM^ as a function of x z - droplet location 
is required.. 

Also the interpolation process of the Exit droplet position (point 4 in 
Figure 4.6-2 ) requires the solution of two algebraic equations viz: 

- droplet trajectory equation; and 

- grid cell face plane equation. 

4 . 7 Interphase Transfer Source Terms 

In the philosophy of the Eulerian-Lagrangean approach, the droplets are regarded 
as source of mass, momentum and energy to the conveying gaseous phase. The source 
terms are incorporated into the gas flow equations, providing the influence of 
the droplet spray on the gas velocity and temperature fields. 

In each grid cell crossed by the droplet trajectory , appropriate interphase 
sources are calculated based on the droplet mass, momentum and energy 
difference between the inlet and the outlet from the cell. 
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The number of particles per unit time, which enter at port i and have an 
initial mass M dj . is given by: 


M dj,- ■ s j h 


(4.6-7) 


where ift d is the total droplet mass inflow rate, X- is the fraction of droplet 
mass which enters at port j, and Y. is the fraction of droplet mass with 
initial diameter D. . The number flow rate of spherical droplets with initial 
diameter D.. , along a given trajectory is calculated as: 


fl = 

y 


6 * d 

u P d D, 3 


(4.6-8) 


where is the droplet density,, 


This value is constant along a droplet trajectory, provided no droplet coalescence 
or shattering takes place. 


Assuming the droplets are spherical, the continuity source term, , representing 
the net efflux rate of droplet mass to the gas phase, is given by: 


M d = 6 ? } : rt ij ^ p d D ^EX " ( p d D ^ INL ^ 

\J 


(4.6-9) 


where the summations are performed over all trajectories crossing the grid cell. 
The "EX" and "INL" subscripts refer to the droplet exit from and inlet into the 
control volume (Figure 4.6-2). The source terms for the remaining dependent 


variables <j> (u, v, w, k, e, n, m. ...) are calculated from the following 

J 

expression: 


4nt,<j> “ Wd - 


(4.6-10) 


where $ d and <Jj represent average <j>- property of the liquid and gas phase in the 
control volume, respectively. Note that the above formula can be conveniently 
linearized as: 


SU INT,<f> Am d • ^d 
SP INT,<f> = " Am d 


(4.6-11) 
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With the assumption that the mass is transferred from liquid to the gaseous 
phase ( rn d > 0). A simple way to calculate droplet average property $ is to 
take an average from the inlet 4^.^ and exit 4> djEX to the control volume and 
then weighted on all droplet sizes. Some of the droplet properties remain 

constant during the integration period and average $ d can be specified directly 
These include: 


^d 1 > for f and m,. 

fu 

^d = C Pd * T sat + H FU » for h 

^d = 0 > for k, e and m CQ 

where Cp d and H py are droplet specific heat, and droplet mass heat of 
combustion. 

4.7.1 Interphase Momentum Transfer 

There are two mechanisms of the momentum transfer between the liquid and vapor 

1) momentum transfer with the associated mass: 

Arfi (u d - u) 


2) frictional momentum transfer: 
f D (u "d 2 " u " 2) 


(4.6-12) 

(4.6-13) 


The second mechanism has not been incorporated in the present version of the 
code and is planned to be included in the future calculations. Nonlinear 
characters of the friction terms requires special linearization practice for the 
momentum equation source terms. The linearization practice for these terms 
should be implemented in the following manner: 

a) the source term can be expressed as: 


S = Su* + — 
u au 


(u - u*) = 


f D ( u d 


u* 2 ) - 2 f 


D 


(u - u*) 


(4.6-14) 


b) appropriate source terms for the u-momentum equation should 
be added as follows: 

SU u = SU u + f D ((I d 2 + 2 I 5 *I 5 * " a * 2 > 


(4.6-15 ) 


SP u " SP u ' 2 f D 
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5. 


BOUNDARY CONDITIONS AND THEIR TREATMENT 


5.1 Introduction 

In order to completely specify the mathematical problem it is necessary to 
supply the conditions at the boundaries of the solution domain for all the 
dependent variables. These conditions are usually either specification of 
the value of the dependent variable at the boundary, or the value of the 
associated flux or a relation between the two. 

Figure 5.1 presents an example of the computational domain and appropriate 
boundary conditions which may in general include: 

1. inlet of gaseous or liquid species (fuel, air, steam, etc.); 

2. exits plane; 

3. center-lines and/or symmetry planes; 

4. periodic boundary conditions; and 

5. solid walls. 




A-A 


Wall 



Figure 5.1. Computational Domain and Boundary Condition Specification 
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For each dependent variable, at all boundaries , appropriate modifications of the 
finite difference equation at the "near boundary" nodes is required. For most 
of the dependent variables the boundary conditions are implemented by modifying 
the source terms SU^ and SP^ of the finite difference equation: 


= 


*d + SU <j> 


Z a 


SP, 


(5.1-1) 


These modification practices are discussed in the following sections. 


5 . 2 Inlet With Specified Flow Rate 

At the inlet the amount of incoming mass flow rate rfij^ through the boundary 
cell face and the incoming <f>g property should be specified. In this case 
the finite-difference coefficient connecting the boundary node to its 
neighboring internal node is set to zero and then the SU^ and SP^ source terms 
are modified as: 



SU cj> + ^IN 


(5.2-1) 



SP , - ih 
$ 


IN 


(5.2-2) 


The pressure correction equation is not modified in this manner. A link 
coefficient with the boundary node is set to zero for p' -equations, 

5 . 3 Symmetry Plane or Axis 

At symmetry plane a zero mass flux ifi^ = 0 is assumed. The modification is 
implemented by setting the appropriate coefficient to zero (Figure 5.2) and 
no modification is required to the source terms. 



Figure 5.2. a $ = 0 at the Symmetry Plane 
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Gas phase injectors, or liquid spray injections with instantaneous evaporation 
assumption, are treated in a similar way. The liquid spray injection for 
which the Lagrangean droplet tracking is employed do not require specification 
boundary conditions in Eulerian meaning. Instead an appropriate initial droplet 
condition should be specified. 


5.4 Exit Boundary 

Where the fluid flows out of the calculation domain, information about most of the 
dependent variables is often not available. However, since it is the process 
occurring in the calculation domain that decide the values of the variables 
which the outgoing fluid will carry, information is not strictly required at 
such boundaries. To treat these boundaries, the boundary coefficients are 
simply set to zero. 

If a fixed exit pressure boundary condition is specified, the velocity normal to 
the exit plane is solved for at the exit boundary and the pressure correction 
p 1 - link coefficient with the exit boundary is calculated. A specified 
pressure correction p' E = 0 at the exit boundary is employed,. 


5 . 5 Periodic Boundary Conditions 

The periodic (cyclic) boundary conditions appear in the circumferential 
direction if the two ends of the calculation domain in the z-direction join 
up with one another,, This can occur in a polar-coordinate direction in which 
the whole angular extent from 0 to 360° is to be considered (Figure 5.3a), or 
when "repetition" is present in the flow pattern in the angular coordinate 
direction (Figure 5.3b)„ 

The general rule is that whenever identical conditions are to be expected at 
z = 0 and z = last z, and finite flow is to be expected at that surface, then 
the boundaries are cyclic. 


In this circumstance, 
c() LB = 


the boundary conditions can be specified as follows: 

■.z 


<f> 


HB 


<f> 


1 


and 


C z = 
L 1 


"HB 


(5.5-1) 


where indices HB and LB denote High Boundary (k = N) and Low Boundary (k = 1), 
C is the convective flux in the z-direction. 
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CYCLIC CONDITIONS FOR d = 120° SECTOR 



Figure 5.3. Grid Notation for Periodic Boundary Conditions in Z-Direction 


5.6 Wall Boundaries 


At the solid walls, the velocity normal to the wall and appropriate convective and 
diffusive fluxes are set to zero. The boundary conditions which can be easily 
specified at the solid walls (Figure 5.4) include: 


v B - 0, 


8(f) 

3y 


= o 

B 


(({> = f, m F|J and m CQ ) 


(5.6-1) 

(5.6-2) 



Figure 5.4. Grid Cell Adjacent to the Solid Wall 




The boundary conditions for velocity components parallel to the wall; for k» 

f\, 

e and h require special "near wall" treatment. 


There are two important features that distinguish near-wall regions from other 
parts of the flow field. Firstly, there are steep gradients of most of the 
flow properties, and secondly, the turbulent Reynolds number is low so that 
the effects of molecular viscosity can influence the shear stresses, 
production dissipation and transport of turbulence energy. 

The vigorous incorporation of these effects requires a prohibitively fine 
finite-difference grid in the vicinity of the wall. 


An alternative approach, however, is available which bridges the near wall 
region and the outer edge of the viscous sublayer by using the "wall functions". 
These are described belowo 


5 . 6 o 1 Wall Functions 

The wall functions described by Patankar and Spalding (Reference 55), Launder 
and Spalding (References 1, 56) and more recently by Launder (Reference 57) 
are derived from experimental and analytical knowledge of the one-dimensional 


Couette flow which exists near the wall. A semi-empirical universal 
of nondimensional distance normal to the wall y + , is: 

+ _ p 6y • u T 


function 


(5.6-3) 


In the above definition 6y is the distance normal to the wall (Figure 5.4 ) 
and u is the "friction velocity" given by: 
t 1/2 

u t = <r> (5 - M) 

In the internal sublayer (y + > 11.63) the velocity variation may be described 
by a logarthmic relationship (see Schlichting (Reference 58)) i.e.: 

u = - a -- An (Ey + ) (5.6-5) 

Is 


where E = 9„793 and k = 0.4187 are experimental ly determined constantSo 
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In both the viscous (y + <^ 11.63) and internal (y + > 11.63) sublayers, the 
shear stress is calculated from the product of effective viscosity y eff - and 
normal velocity gradient 9u/9y, i„e.: 


T w ^eff 9y 
where y 0ff = 


y 

^turb 


for y + £ 11.63 
for y + > 11.63 


(5.6-6) 

(5.6-7) 


Near the wall, the transport equation for the turbulent kinetic energy, k, 
reduces to a balance between the local production and dissipation of k 
(References 1, 56) to give: 


,9iu 2 

y t (jy) = P £: 


(5.6-8) 


The velocity gradient may be replaced from equation (5.6-6) and the dissipation 
rate from: 



(5.6-9) 


to give: 

,, 1/2 . 2 

x = t p k = p u 
w y t 


Hence, it follows from equation (5.6-5) 

r 1/4 . 1/2 - 

p C.. k ' u 


L w 


r 


y 


£n (Ey + ) 


(5.6-10) 


(5.6-11) 


5 „ 6 o 2 Velocity Boundary Conditions 

Equation (5.6-11) is introduced into the finite-difference equation (5.1-1) 
by setting the value of the link coefficient a d (connecting point p with the 
wall node B) to zero, and adding to SP u the term t w . A wall where A wall is the 
cell wall area over which t w acts. For the velocity normal to the wall the 
same process is applied, but the normal shear stress is set to zero. 
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5.6.3 Turbulence Variables 


Due to the steep velocity gradients near the solid walls the assumption of linear 
variation of u is inaccurate and can cause incorrect evaluation of the turbulence 
generation rate G. To overcome this, the generation term near the wall is written as 


r _ 3u 
b k T w 3y 


(5.6-12) 


where t is evaluated according to equation (5.6-11). This is incorporated 
into the difference equation (5.1-1) by setting the link coefficient a d to 
zero and modifying the source terms as follows: 


SU. 


SP, = 


l- 

w 3y 
P 2 C Q k 


w 


3u_ 

9y 


(5.6-13) 


(5.6-14) 


The diffusion of the dissipation rate of turbulence at the wall is a little 
difficult to express. Instead of attempting to calculate r ^ for e, use is 
made of the fact that the length scale £ varies linearly with distance in the 
neighborhood of a wall . Thus, the practice is to "fix" the value of e at the 
near-wall grid point in accordance with: 

e = C D 3/4 k 3/2 (kS) (5.6-15) 


The fixing of e (and similarly for other quantities) can be done by the use 
of the following expressions forSU £ and SP £ of the near-wall point: 


SU = 10 10 * 

e e 

and SP„ = - 10 10 

£ 


(5.6-16) 


When such large values of SU and SP are present, the actual value of a normal 
r wall 1S i mma1:ena l « 


*Instead of 10 3 , any suitable large number may be used, as long as it is 
ensured that the other terms in the finite-difference equation are negligible 
compared to these two terms. 
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5.6.4 Stagnation Enthalpy 

For turbulent flow, expressions have been developed, similar in form, to 
equation 5.6-11, in order to evaluate the heat flux $ across the wall boundary 
layer; the one employed here is due to Jayatillaka (Reference 57) and writes: 


_ 1LX 


{ - Jin 
K 


(Ey + ) 


+ f V 


- 1 d Lw 

dy 


(5.6-17) 


where is the normal temperature gradient. The term, P^, expresses the 
contribution of the laminar sublayer to the total resistance and is calculated 


as: 
P = 


® Q ® 0 

9 - - 1 } {—} 

a h °h 


-1/4 


(5.6-18) 


where a and o are laminar and turbulent Prandtl numbers, respectively. 

v* f I 


For laminar flow (y + < 11.63) the corresponding expression for the heat flux 
is: 


! T w 

°o dy 


(5.6-19) 


Equation (5.1-1) for h is modified by breaking the link between the near-wall 

nodes and adding to SU h the term $ w . A ^ , where A is the cell-wall area 

through which is transferred, 
w 


5 . 7 Boundary Conditions for the Radiation Equations 

A provision of zero gradient of net radiative heat fluxes at the calculation 
domain boundaries is provided as the default case. This is achieved by 
making relevant boundary link coefficients to zero. Other boundary specifications 
can be handled by the modifications of the near-boundary source terms. This 
is outlined below. 
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Making the boundary coefficient zero implies that ( Td R /dx ) for the boundary 
cell surface is made zero. This can be seen to be equal to — Q /2 from equation 
3.6-3. 


Thus, the required source-term modification for the cell close to the boundary 
can be stated as the inclusion of - Q /2 for the boundary as an additional source 

X 

term. For various particular cases, this will involve the following additional 
contributions to SU and SP. 

Symmetry Plane 

At a symmetry plane, Q is zero by definition. Hence no modification of SU 
and SP is needed. 

Non-Reflecting Boundary 

If the outgoing radiation leaves the calculation domain without reflection 
and if the incoming radiation, equal to, say L, is given, from the definitions 
of R x and Q x , 

-Q x /2 = L - R x (5.7-1) 

Thus the additional SU should be L (which is given) and the additional SP 
should be -1. 

Wall Boundary 

If is the emissivity of the wall and the black-body emissive power at 

the wall temperature, the flux leaving the wall, say L, is given by: 

L = e w E w + 0 - Si> K (5 - 7 " 2 > 

emitted reflected 


Again, via the definitions of Q and R , the following relation is obtained 

X X 

from equation (5.7-2): 


'2 - e, 


-) E„ - ( 


2 - e, 


-) K 


(5.7-3) 
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Thus, {e w E w /(2 - e w )} becomes the additional SU while {-e w /(2 - e w )} is the 
additional SP. 

The y- and z--direction fluxes are treated in a similar manner. It should, 
however, be remarked that for the radial direction the additional source 
term for the y-direction flux is (-rQ^/2) , where r is the radius of the 
boundary surface. 
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NOMENCLATURE 


\|»\» • * * 


a E’ a W’ ** * 


ox 

a j 

B 

b. 

J 

C 

C 1 * C 2* C D 
Cn 


c j 


DU.DV.DW 




f 

G k 

g 

h j 

% 

i 


j 

K 


Area of the grid cell (W-face, E-face, ...) 

Absorption coefficient for radiation 

Finite difference link coefficient (E-east, W-west, N-north, .„.) 

Oxygen mass fraction in oxidizer 

Polynomian coefficient for cp-calculation 

Eddy-Break-Up model constant 

Polynomian coefficient for cp-calculation 

Convective flux 

Turbulence model constants 

Specific heat of mixture at constant pressure 

Polynomian coefficient for cp-calculation 

Droplet diameter 

Pressure-di fference coefficients 

Polynomian coefficient for cp-calculation 

Activation energy in the Arrhenius reaction rate law (Chapter 3.7); or 
Black body emissive power (Chapter 3.6); or 
A constant in the law of the wall (Chapter 5) 

Mixture fraction 

Generation rate of turbulence energy 
Concentration fluctuation 
Heat of combustion for j-th species 
Stagnation enthalpy 

False inertia term (Chapter 4.5); or 

Radiation flux in the positive x-direction (Chapter 3.6) 

Radiation flux in the negative x-direction 
Radiation flux in the positive r-direction 
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k 

L 

M 


m 

N 

NS 

n 

P 


Qx’Vz 


R J 


R x ,R y ,R z 


r 

S 

SU, SP 
s 


T 

t 


-■ Kinetic energy of turbulence 

Radiation flux in the negative r-direction 
~ Mixture molecular weight (Chapter 3.3); or 
Radiation flux in the positive z-direction (Chapter 3.6); or 
Mass of the control volume (Chapter 4.3) 

Molecular weight of species j 
- Mass fraction of species j 

- Coefficient in hydrocarbon composition (C n H ) specification 

- Radiation flux in the negative z-direction 

- Number of species participating in the mixture composition 

- Coefficient in hydrocarbon composition (C H m ) specification 

- Resistance of the laminar sublayer (Chapter 5.6); or 

Pre-exponential factor in Arrhenius reaction rate expression 
(Chapter 3.7) 

- Pressure 

- Pressure correction 

- Net radiative heat fluxes in the x, y and z (or 0) directions 

- Mass rate of creation of species j by chemical reaction 

- Universal gas constant 

- Composite radiation fluxes (dependent variables) in the x, y, 
z (or 0) directions 

- Distance from axis of symmetry 

- Source term 

- Parts of linearized source term 

- Mass ratio of stoichiometric osidant/fuel proportions (Chapter 3.7); or 

- Radiation scattering coefficient (Chapter 3.6) 

- Absolute temperature 

- Time 
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U, V, w 
U 

W 

X, y, Z 

+ 

y 


- Velocity components in the x, y and z (or 0) directions 

- Velocity vector 

- Angular velocity (W = rw) 

- Coordinate distances 

- Dimensionless distance from the wall 


Greek Symbols 

a - Relaxation factor 

^ = m fu - ”ox /s 

r - Diffusion coefficient 

Ax,Ay,Az - x-direction, y-direction and 6 (or z) direction lengths of a 
control volume 

Sx,Sy,6z - x-, y- and z-direction distances between the node points 
e - The dissipation rate of turbulence 

- Wall emissivity 

0 - Coordinate distance 

k - Von Karman constant 

U - Viscosity 

p - Density 

a - Laminar Prandtl/Schmidt number (Chapter 3.8); or 

Stefan-Bol zman constant (Chapter 3.6) 
c p - The general dependent variable 

- Spherical angle 


Subscripts 

A - Air inlet 

d - Droplet 

eff - Effective value 

E - East side (x + ) neighbor 

EBU - Eddy-break-up 

fu - Fuel 
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F 

- Fuel inlet 

H 

- High size (z + ) neighbor 

ft 

- Stagnation enthalpy 

int 

- Interphase value 

j 

- Species j 

k 

- Kinetic energy of turbulence 

L 

- Low size (z") neighbor 

1 

- liquid 

l 

- Laminar 

N 

- North side (y + ) neighbor 

ox 

- Oxygen 

pr 

- Products 

stoich 

- Stoichimetric 

S 

- South side (y~) neighbor 

t 

- Turbulent 

V 

- Vapour 

W 

- West side (x - ) neighbor 
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